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Abstract 

The goal of this paper is to generalize the well-balanced approach for non-equilibrium 
flow studied by Wang et al. [26] to a class of low dissipative high order shock-capturing filter 
schemes and to explore more advantages of well-balanced schemes in reacting flows. The 
class of filter schemes developed by Yee et al. [30], Sjogreen & Yee [24] and Yee & Sjogreen 
[35] consist of two steps, a full time step of spatially high order non-dissipative base scheme 
and an adaptive nonlinear filter containing shock-capturing dissipation. A good property 
of the filter scheme is that the base scheme and the filter are stand alone modules in design- 
ing. Therefore, the idea of designing a well-balanced filter scheme is straightforward, i.e. , 
choosing a well-balanced base scheme with a well-balanced filter (both with high order). 
A typical class of these schemes shown in this paper is the high order central difference 
schemes/predictor-corrector (PC) schemes with a high order well-balanced WENO filter. 
The new filter scheme with the well-balanced property will gather the features of both 
filter methods and well-balanced properties: it can preserve certain steady state solutions 
exactly; it is able to capture small perturbations, e.g., turbulence fluctuations; it adap- 
tively controls numerical dissipation. Thus it shows high accuracy, efficiency and stability 
in shock/turbulence interactions. Numerical examples containing ID and 2D smooth prob- 
lems, ID stationary contact discontinuity problem and ID turbulence/shock interactions 
are included to verify the improved accuracy, in addition to the well-balanced behavior. 

Key words: High order filter methods, WENO schemes, well-balanced schemes, non- 
equilibrium flow, chemical reactions, ID turbulence/shock interactions. 


1 Introduction 

Recent progress in the development of a class of low dissipative high order filter schemes for 
multiscale Navier-Stokes and magnetohydrodynamics (MHD) systems [30, 37, 24, 32, 23, 33, 
34, 35] shows good performance in multiscale shock/turbulence simulations. 

The highly parallelizable high order filter methods consist of two steps, a full time step 
of spatially high order non-dissipative (or very low dissipative) base scheme and an adaptive 
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multistep filter. The nonlinear filter consists of the products of wavelet based flow sensors and 
the dissipative portion of high order shock-capturing schemes. The built-in flow sensors in 
the post processing filter step can control the amount and types of numerical dissipation only 
where needed, and leave the rest of the flow region free from numerical dissipation. Only the 
filter step might involve the use of flux limiters and approximate Riemann solvers as stabilizing 
mechanisms to remove Gibb’s phenomena related spurious oscillations resulting from the base 
scheme step at locations where needed. The more scales that are resolved by the base scheme, 
the less the filter is utilized, thereby gaining accuracy and computational time. The adaptive 
numerical dissipation control idea is very general and can be used in conjunction with spectral, 
spectral element, finite element, discontinuous Galerkin, finite volume and finite difference 
spatial base schemes. The type of shock-capturing schemes used as nonlinear dissipation can be 
the dissipative portion of any high resolution TVD, MUSCL, ENO, or WENO shock-capturing 
methods [30, 11, 21]. By design, the flow sensors, spatial base schemes and linear and nonlinear 
dissipation models are stand alone modules. Therefore, a whole class of low dissipative high 
order schemes can be derived at ease. 

In the recent paper by Wang et al. [26], well-balanced finite difference WENO schemes 
and second-order TVD schemes were studied for chemical non-equilibrium flows, extending the 
well-balanced finite difference WENO schemes for shallow water equations in [27, 28]. A well- 
balanced scheme [13], which can preserve certain nontrivial steady state solutions exactly, may 
help minimize some of the oscillations around steady states. It was also shown in [26] that 
the well-balanced schemes capture small perturbations of the steady state solutions with high 
accuracy. While general schemes can only resolve perturbations at the level of truncation error 
with the specific grid, well-balanced schemes can resolve much smaller perturbations, usually 
at 1% or lower of the main steady state flow. 

In this paper the same approach will be applied to construct a high order well-balanced filter 
scheme for one temperature non-equilibrium flow with reaction terms. The multi-dimensions 
hyperbolic system of conservation laws with source terms (also called a balance law) 

U t + V-F(U) = S(U) (1) 

is considered, where U is the solution vector, F(U) is the convective flux and S(U) is the 
source term. For this type of flow the space variable x does not appear explicitly in the source 
term. Thus, the construction of well-balanced schemes can easily go from one-dimension to 
multi-dimensions. 

The designing of well-balanced filter schemes is to choose a well-balanced base scheme and 
a well-balanced filter part. Then, the filter scheme is almost well-balanced everywhere except 
at the interfaces of the filtered region and the non-filtered region (see Sec. 5). Note that in 
this paper, a ‘well-balanced filter scheme’ refers to such almost well-balancedness. For the filter 
scheme without the flow sensor, the resulting filter scheme is well-balanced. 

The choice of the sensor will not destroy the well-balanced properties. It has been shown 
in the previous work [26] that linear schemes, second-order Predictor-Corrector (PC) scheme 
[31, 14] with TVD filters (such as Harten-Yee TVD filter [29, 30]) and WENO-Roe schemes are 
well-balanced for certain steady state solutions with zero velocity. A well-balanced WENO- 
LF scheme has also been constructed for this type of steady state solutions. High order 
PC schemes are linear schemes and thus well-balanced. Therefore, the new filter schemes 
CENTVDfi/CENWENOfi or PCTVDfi/PCWENOfi which utilize central/PC schemes as the 
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base schemes and the Harten-Yee TVD filter or the balanced WENO schemes as the filter will 
be well-balanced. 

In this paper, only the “zero velocity” steady state of the reacting flow will be considered in 
the numerical tests. Zero velocity steady state means the velocity is zero which will also imply 
that the flow has constant pressure and stay in chemical equilibrium. It will be shown that 
same as the well-balanced WENO schemes, well-balanced filter schemes give machine round-off 
errors regardless of the mesh sizes for the stationary steady state solutions of the reactive flow. 
Consequently, they can resolve small perturbations of such steady state solutions well with very 
coarse meshes. 

Since the regular high order low dissipative filter schemes are designed for slrock/turbulence 
interactions and the well-balanced schemes can capture small perturbations of the steady state 
solutions with high accuracy, the new well-balanced filter schemes take the advantages of both, 
thereby making them well suited for computations of turbulent fluctuations on a mainly steady 
flow field. 

The outline of the paper is at follows: in Sec. 2, the governing equations and physical 
model are described. The high order filter schemes and well-balanced schemes are reviewed in 
Sec. 3 and 4. Then the proposed high order well-balanced filter scheme is introduced in Sec. 
5. Numerical examples will be shown in Sec. 6. Finally, Sec. 7 is the conclusion and future 
work. A brief description of high order PC schemes and the considered time discretization are 
presented in the Appendix. 


2 Governing equations and physical model 

Considering an atmospheric entry flow in chemical nonequilibrium and thermal equilibrium, 
the thermodynamic properties account for excitation of the electronic states for the atoms and 
molecules, using the rigid- rotor harmonic-oscillator approximation for the molecules [15, 18]. 

Assuming neither dissipative effects nor radiation, the considered physical model is a system 
of hyperbolic conservation laws with source terms denoted by 

U t + V • F{U) = S(U), (2) 
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where n s is the number of species, p s , the mass density of species s, u, the velocity vector, and 
e, the internal energy per unit mass of the mixture. The mixture mass density is defined as 
p = P s > an< 4 the pressure p is given by the perfect gas law 
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where R is the universal gas constant, and M s , the molar mass of species s. The temperature 
T can be found from the total energy 


P e = ^2pse s (T) + -p|u| : 

S— 1 


( 5 ) 


where the internal energy e s of species s is a function of temperature 


es(T) = ei(T) + ef(T) + e 
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where H a is the set of atoms and H p is the set of molecules. The translational energy of species 
s is given by 

e T (T) = -—T 
A ’ 2 M s ’ 

and the electronic energy contribution by 
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where quantities gf n and 9f n stand for the degeneracy and characteristic temperature of the 
electronic level n of species s. The number of electronic levels retained is limited for mathe- 
matical and physical standpoints. The partition function leading to thermodynamic properties 
diverges when all levels are accounted for. The maximum number of electronic levels of each 
atom and molecule is progressively increased up to a correspondence between the values of 
computed enthalpies and accurate reference tables. For molecule s, the rotational energy is 
assumed to be described by means of a rigid rotor model 


l (T) = 


RT 
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and the vibrational energy, by means of a harmonic oscillator model 
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where the quantity 6Y stands for the vibrational characteristic temperature of the diatomic 
molecule. To account for the energy released in the gas by chemical reactions between the 
species, a common level from which all the energies are measured is established by using the 
formation enthalpy ef at 0°K. 

The source term S(U) describes the chemical reactions occurring in gas flows which result 
in changes in the amount of mass of each chemical species. In the general case there are J 
reactions of the form 
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where v' x • and i/" • are respectively the stoichiometric coefficients for the reactants and products 
of species i in the jth reaction. For non-equilibrium chemistry the rate of production of species 
i due to chemical reaction may be written as 




i,3 


3 = 1 




M, 


i = 1, 


(13) 


For each reaction j the forward and backward reaction rate coefficients, kfj and kbj are assumed 
to be known functions of temperature. The forward reaction rate coefficient is given by an 
Arrhenius law. Following microreversibity the backward rate coefficient is obtained from the 
expression kfj = kbj/K e j, where the equilibrium constant for the jth reaction is given by the 
relation 

InKej = ~ V 's,j) m s9s{Pref,T)\, (14) 

kb t “ 

S= 1 

where the reference pressure p re f = 1 Pa. The Gibbs free energy g s of species s is a function 
of pressure and temperature 

g s (p,T) = gJ(p,T) + g?(T), s G H a , (15) 

9 s(p,T) = gJ(p,T)+g?(T)+g?(T)+gY(T), s G H p . (16) 


The translational Gibbs free energy is obtained from 
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where the symbol hp stands for Planck’s constant, and Na for Avogadro’s number. The 
electronic Gibbs free energy reads 


gf(T) = --T7- in 
as v ) M 
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For the diatomic molecule s, the rotational Gibbs free energy is 
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where symbol a s stands for the steric factor. The vibrational Gibbs free energy is obtained 
from the relation 

1 — exp 
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3 Description of high order filter methods 

For simplicity, the numerical methods for the ID case is described. Denote A = dF/dU , the 
Jacobian matrix of the flux in Eq. (3). The eigenvalues of A are 

(a 1 , . . . , a m ) = (u, . . . ,u,u + c,u — c), (21) 
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where m is the number of components of vector U, m = n s +2 in ID case, c is the frozen speed of 
sound defined by the expression c 2 = (n+l)p/p with quantity n = (X^”=i PsR/M s )/(Y^Li PsC v ,s) 
based on the species specific heat c„ iS = de s /dT . Denote R as the matrix whose columns are 
eigenvectors of A (not to be confused with the R in Eq. (4)). Let a l . +1 , 2 , Rj+\/2 denote the 
quantities a 1 and R evaluated at some symmetric average of Uj and Uj+i, such as Roe’s average 
[22], Define 

a j+i/2 = Rjh/ 2 (Uj+i -Uj) (22) 

as the difference of the local characteristic variables in the x direction. 

The low dissipative high order filter scheme developed by Yee et al. [30] suggests using the 
artificial compression method of Harten [7] as a flow sensor to limit the amount of numerical 
dissipation that is inherent in a scheme. Subsequently, Sjogreen and Yee [24], Yee and Sjogreen 
[33, 35] introduced a wavelet decomposition of the data to determine the location where numer- 
ical dissipation is needed. The considered filter method contains two steps, a high order low 
dissipative spatial base scheme step (not involving the use of approximate Riemann solvers or 
flux limiters) and a multistep filter (usually involving the use of approximate Riemann solvers 
and flux limiters). The nonlinear filter consists of the product of a wavelet (WAV) sensor [24] 
and the nonlinear dissipative portion of a high-resolution shock-capturing scheme. 

We will briefly review the high order filter schemes in this section. For more details, we 
refer the readers to [30, 24, 33, 35]. 


3.1 High order spatial scheme step 

The first step of the numerical method consists of a time step via a high order non-dissipative 
spatial and high order temporal base scheme operator L* . After the completion of a full time 
step of the base scheme, the solution is denoted by U* 

U* = L*(U n ), (23) 


where U n is the numerical solution vector at time level n. 

The high order non-dissipative spatial base schemes could be the standard central schemes, 
the centered compact scheme or the predictor-corrector (PC) schemes [31, 14]. 

For strong shock interactions and/or steep gradient flows, a small amount of high order 
linear dissipation can be added to the base scheme step to reduce the time step constraint and 
stability. For example, an eighth-order linear dissipation with the sixth-order central scheme 
to approximate F(U) X is written as 


dF 

dx 


Doe F 3 + d(Aa0 7 (D + D„) 4 [/„ 


(24) 


where Doe is the standard sixth-order accurate centered difference operator, and D + D~ is the 
standard second-order accurate centered approximation of the second derivative. The small 
parameter d is a scaled value (e.g., spectral radius of A(U)) in the range of 0.00001-0.0005, 
depending on the flow problem, and has the sign which gives dissipation in the forward time 
direction. See Appendix for more choices of spatial base schemes. 
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3.2 Adaptive nonlinear filter step (discontinuities and high gradient 
capturing) 

After the completion of a full time step of the high order base scheme, the second step is to 
adaptively filter the solution by the product of a “wavelet sensor” and the “nonlinear dissipative 
portion of a high-resolution shock-capturing scheme” (involving the use of flux limiters). The 
final update step after, e.g., the nonlinear filter step can be written as 

v r = u s-^ x Kv - ■ < 25 > 

Here, the filter numerical fluxes 

Hj+ 1/2 = Rj+i/ 2 Hj+i/ 2 ■ (26) 

Denote the elements of the vector Hj+ 1/2 by /i( +1 , 2 , l = 1, 2, . . . m. The nonlinear portion of 
the filter liA +1 , 2 has the form 

^i/2 = \{s N ) l j+1/ 2^ l j + y2)- (27) 

Here, (s JV )(. +1 y 2 is the sensor to activate the higher order nonlinear numerical dissipation <?r +1 / 2 - 
(s N ) l j +1 / 2 is designed to be zero or near zero in regions of smooth flow and near one in regions 
with discontinuities. (■s Ar )) +1 / 2 var i es from one grid point to another and is obtained from 
a wavelet analysis of the flow solution [24]. The wavelet sensor can be obtained from the 
characteristic variables for each wave or a single sensor for all waves, based on pressure and 
density. The latter one is used in the paper. 

The dissipative portion of the nonlinear filter <j>j +1 / 2 = 9j+i / 2 — ^j+ 1/2 the dissipative 
portion of a high order high-resolution shock-capturing scheme for the local 1-th characteristic 
wave. Here t/) +1 / 2 anc ^ ^j+ 1/2 are numerical fluxes of the uniformly high order high-resolution 
scheme and a high order central scheme for the 1-th characteristic wave, respectively. 

For the numerical examples, two forms of nonlinear dissipation </>) +1 ^ 2 were considered, 
namely: 

• Dissipative portion of balanced WENO schemes (see Sec. 4). It is obtained by taking 
the full WENO scheme and subtracting the central scheme, such as, WEN05-Do6 and 
WEN07-A)8. 

• Dissipative portion of the Harten-Yee TVD scheme [29, 30]. 

Remark 1. In [35], Yee and Sjogreen proposed and studied both linear and nonlinear filters, 
where the linear filters refer to the standard spectral filter, compact filter and non-compact high 
order linear numerical dissipation. In our paper, only nonlinear filters, especially the dissipative 
portion of WENO and TVD schemes are considered. 

Remark 2. Yee and Sjogreen also did comparisons of applying the filters between “after each 
Runge-Kutta stage” and “after a full time step”. Their research indicated that there is no 
advantage of applying the filters “after each Runge-Kutta stage”. In addition, “after a full 
time step ” is extremely efficient since only one Riemann solve per time step per dimension is 
required. 
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3.3 Flow sensor by multiresolution wavelet analysis of the computed 
flow data 

A general description of how to obtain different flow sensors (e.g., (/)- +1/2 ) by multiresolution 
wavelet analysis of the computed flow data can be found in Sjogreen and Yee [24] and Yee and 
Sjogreen [32]. 

The wavelet flow sensor estimates the Lipschitz exponent of a grid function fj (e.g., the 
density and pressure). The Lipschitz exponent at a point x is defined as the largest a satisfying 


sup 

h^O 


I f(x + h) - f{x) | 
h a 


<C, 


(28) 


and this gives information about the regularity of the function /, where small a means poor 
regularity. For a C 1 wavelet function w with compact support a can be estimated from the 
wavelet coefficients, defined as 

Wm,j = (f,1pm,j) = j f(x)lpm,j(x)dx, (29) 

where 

Vw = 2 ro V>(^=r) (30) 

is the wavelet function tp m ,j on scale m located at the point j in space. In practical computations 
we have a smallest scale determined by the grid size. We evaluate w m j on this scale, too, and 
a few coarser scales, too + 1, too + 2, and estimate the Lipschitz exponent at the point jo by a 
least square fit to the line [24] 


max log 2 \w m j | = ma + c. (31) 

j near j 0 


A continuous function f(x) has a Lipschitz exponent a > 0. A bounded discontinuity (shock) 
has a = 0, and a Dirac function (local oscillation) has a = — 1. Large values of k can be used 
in turbulent flow so that large vortices or vortex sheets can be detected. 

For more details about the wavelet and the flow sensor, we refer the readers to [24, 35] . 


4 Description of well-balanced methods 

A well-balanced scheme refers to a scheme that preserves exactly specific steady state solutions 
of the governing equations. In the previous work [26] the linear schemes, WENO-Roe scheme, 
Harten-Yee TVD and the Predictor-Corrector TVD schemes (with zero entropy correction) are 
proven theoretically and numerically to be well-balanced schemes for the non-equilibrium flow 
Eq. (2) with zero velocity steady states. We will briefly review the idea of the well-balancedness 
in this section. 

For the general 1-D system balance law 


U t + F(U,x) x = S(U,x), 


(32) 



the steady state solution U satisfies 

f\U,x) x = s l (U,x), l = l,...,m, (33) 

where f l and s l are the Zth elements of the vectors F(U,x) and S(U,x). 

A linear finite-difference operator D is defined to be one satisfying D(afi +6/2) = aD(fi) + 
bD(f 2 ) for constants a, b and arbitrary grid functions /1 and f 2 . A scheme for Eq. (32) is said 
to be a linear scheme if all the spatial derivatives are approximated by linear finite-difference 
operators. 

As proved in [27], under the following two assumptions regarding Eq. (32) and the steady 
state solution of Eq. (33), linear schemes with certain restrictions are well-balanced schemes. 
Furthermore, high-order nonlinear WENO schemes can be adapted to become well-balanced 
schemes. 

Assumption 1. The considered steady state preserving solution U of Eq. (33) satisfies 

r s (U, x ) = constant , s = 1,2,... (34) 

for a finite number of known functions r s (U, x). 

Assumption 2. Each component of the source term vector S(U,x) can be decomposed as 

s\U,x) = ^2,Ti{r 1 {U,x),r 2 {U,x),...)t\{x), l = (35) 

i 

for a finite number of functions Ti and U, where Ti could be arbitrary functions of r s (U,x), and 
Ti and ti can be different for different s l (U,x). (Here ti is not to be confused with the time “t” 
indicated on all previous conservation laws.) 

Now consider the non-equilibrium flow Eq. (2). First, since S(U, x) = S(U ), all the t'^x) = 1. 
Next, when the flow is in the steady state, the chemistry is in equilibrium and thus the source 
vector S(U) = 0. Therefore, two assumptions are easily satisfied by taking 

r i(U) = s % (U), i= 1, . . . ,m. (36) 

Furthermore, linear schemes (such as central and PC schemes) and WENO-Roe scheme are 
naturally well-balanced for such steady state solutions of Eq. (2) (see [26]). 

A well-balanced finite difference WENO-LF scheme can be constructed with a limiter A in 
the Lax-Friedrichs flux splitting 


f ± ( u ) = ^(/O); ± aXu )• ( 37 ) 

A is close to 0 or 1 according to the solution which is in steady state or away from steady state. 
In particular, A is constructed by 


\ \ ’ \r 1 (U i+ i,Xi + i)-r 1 (Ui,Xi)\‘‘ ! + \r 1 (Ui-i,Xi-i)-r 1 (Ui,'r.'i\'- ! 


-\~£ . 


rin (l, 


(\r2(Ui+i,Xi+i)-r2(Ui,Xi)\ + \r2(Ui-i,Xi-i)-r 2 (Ui,Xi)\) 
r2(Ui +1 ,Xi + i)-r 2 {Ui,Xi)\‘ 2 +\r2{Ui-uXi-i)-r 2 {Ui,Xi)\'' i +£ 


(38) 
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where £ is a small number to avoid zero in the denominator and we take it as 10~ 6 in the 
computations. Near the specific steady state, the differences in ry shown in (38) are close to 
zero. A will be near zero when all these differences are small compared with e. A is near one if 
the solution is far from the steady state, since the differences in ry shown in (38) are now on the 
level of 0( Ax) and much larger than £, and then the scheme is the regular WENO-LF scheme. 
The limiter does not affect the high order accuracy of the scheme in smooth region for general 
solutions of Eq. (2). In the specific steady state, since all the ry are constants, A becomes zero 
and then the scheme maintains the exact solutions for such steady state. 

We remark that the functions ry in the limiter (38) are used to distinguish the states be- 
tween steady and unsteady. They are not necessarily to be the same as in Assumption 1, 
but they must be a necessary condition for the steady states. For example, in the consid- 
ered “zero velocity” steady state, zero velocity will imply zero source terms and constant 
pressure etc. u = 0 is a necessary condition for zero velocity steady state. Thus taking 



(\v(Ui + 1 ,Xi+ 1 )\ + \v(Ui,x i )\) 2 
v(Ui+i,Xi+i)\ 2 +\v(Ui,Xi)\ 2 +e 


j also works in the code implementation. 


The dissipative portion of the TYD schemes has the form 



^+1/2) ^ H+ 1 / 2 ) 2 ] ( 


* 7 + 1/2 



(39) 


where ^ 

v )+ 1/2 = Aa- a i+i/2’ 

a j+i /2 and a j+ 1/2 are defined in (21) and (22). 

The function ip( z ) is an entropy correction to \z\ (see [29]) with 





M M><5i 

(z 2 + Sl)/2Si \z\ < <5r ’ 


(41) 


where (5i is the entropy fix parameter (see [36] for a discussion). Q l j + i / 2 is an unbiased limiter 
function which can be 


Q)+ 1/2 = minmod(a*-_ 1/2 , a*- +1/2 ) + minmod(a' +1/2 , a l j+3/2 ) - a l j+1/2 (42) 

with 

minmod(a, b ) = sgn(a) • max{0, min[|a|, b sgn(a)]}. (43) 

It has been proven in [26] that the zero velocity steady state solution, -Ry+i/ 2 -Hj+i /2 in 
Eq. (26) maintains zero and thus can be used as the filter part for the well-balanced filter 
schemes. 


5 High order well-balanced filter scheme 

The construction of high order well-balanced filter schemes is rather straightforward. The 
first step is to choose any well-balanced scheme, such as central scheme (CENa;), PCcc (here 
x denotes the number 2, 4, 6 or 8) or any linear scheme as base scheme. The second step is 
to choose a well-balanced filter, such as the dissipative portion of the TVD scheme (39) or the 
high order well-balanced WENO scheme. 
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Here, we would like to remark that these constructed filter schemes are well-balanced, except 
at the interface between the filtered and non-filtered regions. Because in the interface of these 
two regions, the numerical fluxes get information from different schemes (base scheme part and 
filter part), the schemes will not be well-balanced at those interface cells. This is not a serious 
concern, since the interface is only a small portion of the whole computational domain. Also, 
since the filter is turned on only at the shock region, the transition region of the shock is usually 
far away from the considered zero velocity steady state. Thus, there is no need to require the 
schemes to be well-balanced at the interfaces. 

Also note that the linear dissipation part d(Ax) 7 (D + D_) 4 Uj in the base scheme (24) cannot 
preserve the steady state solutions. Similar to LF flux, since there are no assumptions on the 
density functions, the dissipation d(Ax) 7 (D + D_) 4 Uj may produce non-zero values in the steady 
states. Here, the same idea of constructing well-balanced WENO-LF schemes is applied, i.e., 
multiplying a limiter A (38) to the linear dissipation part to turn off the linear dissipation in 
the steady state area. Since the linear dissipation is only needed for stability concern before 
reaching steady state, numerical tests show that turning it off by the limiter A does not affect 
the stability of the solution. With the limiter A, the filter schemes will have no linear dissipation 
in the steady state and thus will maintain the exact steady state solutions. 

In this paper the considered well-balanced filter schemes are low order central filter schemes 
with TVD filter (CEN2TVDfi and CEN4TVDfi), and high order central filter schemes with 
balanced WENO filter (CEN6WEN05fi and CEN8WEN07fi). Similar considered PC filter 
schemes are PC2TVDfi, PC4TVDfi, PC6WEN05fi and PC8WEN07fi. For the same order PC 
& central filter schemes, the accuracies are similar. Comparing to central filter schemes, PC 
filter schemes allow a larger CFL number for time integration. However, the PC filter schemes 
only allow 1st & 2nd-order time discretizations. The central filter schemes allow a wider class 
of time discretizations. 


6 Reaction model and test cases 

In this section, the gas model for sub-orbital Earth reentries comprising five species N 2 , O 2 , 
NO, N, and O is described. Then, different numerical tests of the considered high order well- 
balanced filter schemes for one- and two- dimensional reacting flows are performed. The first 
example is to numerically verify that the constructed filter schemes are well-balanced by time- 
marching on a nontrivial steady state. In this test the well-balanced filter schemes will show 
round-off numerical errors for a specific steady state solution. The second example is a small 
perturbation over the steady state. We can observe the well-balanced filter schemes showing 
their advantage in resolving the perturbations in very coarse meshes. For ID numerical tests, 
we show three additional examples involving shocks. The first one is a stationary contact dis- 
continuity problem, where the left and right states of the discontinuity are both in equilibrium. 
We will show that if there are small perturbations on the two sides of the discontinuity, the 
well-balanced schemes can capture them very accurately. The second shock example is a ID 
turbulence/shock interaction problem where only the right state of the shock is in equilibrium. 
If there are small perturbations on the right of the shock, the well-balanced schemes will well 
resolve them, then when the shock passes through those perturbations, well-balanced schemes 
will have more accurate results than the non well-balanced scheme on the left of the shock. The 
third example is a shock tube problem to test the shock-capturing capability of the considered 
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schemes. There numerical test cases are to demonstrate that well-balanced schemes will not 
destroy the non-oscillatory shock resolution away from the steady state. 


6.1 Reaction model 

The air mixture comprises 5 species, N 2 , O 2 , NO, N, and O, with elemental fractions 79% 
for nitrogen and 21% for oxygen. The spectroscopic constants used in the computation of the 


species thermodynamic properties (dY , 0f n , df, gf in Sec. 2) and the formation enthalpies 
are obtained from Gurvich et al. [6]. The chemical mechanism comprises three dissociation 
recombination reactions for molecules 

N 2 + M ^2N + M, 

(see [20]) 

(44) 

0 2 + M ^2 0 + M, 

(see [20]) 

(45) 

NO + M ^ N + O + M, 

(see [19]) 

(46) 

where M is a catalytic particle (any of the species N 2 
reactions for NO formation 

, 0 2) NO, N 

, and O), and two Zelclovich 

N 2 + 0 ^ NO + N, 

(see [2]) 

(47) 

0 2 + N ^ NO + O, 

(see [3]). 

(48) 


6.2 One dimensional numerical results 

6.2.1 Well-balanced test 

The purpose of the first test problem is to numerically verify the well-balanced property of the 
proposed filter schemes. The special zero velocity stationary case with 

T = 1000 x (1 + 0.2 sin(7ra;)) K, p = 10 5 N/m 2 , u = 0 m/s, (49) 

is considered. The initial composition is based on the local thermodynamic equilibrium (LTE) 
assumption. Given Ecp (49) and the source term S(U) = 0, each species is uniquely determined. 

Eq. (49) is chosen as the initial condition which is also the exact steady state solution, and 
the results are obtained by time-accurate time-marching on the steady state. The computational 
domain is [—1,1]. The L 1 relative errors of temperature at t = 0.01 (about 1000 time steps 
for N = 100 grid points) are listed in Table 1. The L\ relative error is measured to be the 
difference between the exact solution Eq. (49) and the numerical solution divided by the L% 
norm of the exact solution. 

Table 1 shows that the considered high order central filter schemes and PC filter schemes 
are well-balanced because they produce errors at the level of machine round-off errors in double 
precision. 
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Table 1: L 1 relative errors for temperature by central/PC filter schemes at t = 0.01. 


N 

error 

error 

error 

error 


CEN2TVDfi 

CEN4TVDfi 

CEN6WEN05fi 

CEN8WEN07fi 

50 

3.84E-11 

3.84E-11 

3.79E-11 

3.67E-11 

100 

3.79E-11 

3.79E-11 

3.68E-11 

3.62E-11 


PC2TVDfi 

PC4TVDfi 

PC6WEN05fi 

PC8WEN07fi 

50 

3.83E-11 

3.76E-11 

3.69E-11 

3.63E-11 

100 

3.88E-11 

3.85E-11 

3.81E-11 

3.78E-11 




Figure 1: Small perturbation of velocity results by filter schemes: e = 10 3 x sin(7rx). Left: 
central filter schemes; Right: PC filter schemes. 


6.2.2 Small perturbation 

The following test problem will demonstrate the advantages of well-balanced schemes through 
the problem of a small perturbation over a stationary state. 

The same stationary solution, Eq. (49), is considered. A small perturbation e = 10’ - 3 x 
sin(7nr) is added to the velocity, i.e., 

u' = u + e. (50) 

The other quantities are kept unperturbed. Figs. 1 show the velocities by central and PC filter 
schemes of orders 2, 4, 6 and 8 at t = 0.1. The reference results are computed by fifth order 
WENO-Roe with 1200 points and are considered to be “exact”. 

The results show that all the considered high order well-balanced filter schemes can capture 
the small perturbation well in a very coarse mesh. Especially for the schemes with order higher 
than 2, only 50 points are used. However, the non well-balanced schemes behave in a very 
oscillatory fashion, such as WENO-LF, without the limiter lambda in (38), with 200 points 
(Fig. 2). They can only resolve the solution when the mesh is refined enough such that the 
truncation error of the scheme is much smaller than the perturbation. 
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Figure 2: Small perturbation of velocity results by WENO-LF scheme: e = 10 3 x sin(7rx). 
WENO-LF 200 points: dash-dot; Reference 1200 points: solid). 


6.2.3 ID stationary contact discontinuity problem 

The third example is a ID stationary contact discontinuity problem on the domain [—5,5]. 
A stationary contact discontinuity is located at x = 0. The flow contains zero velocity and 
constant pressure 20N/m 2 everywhere. The temperature has an initial condition 

_ f 500 x (1 + 0.1 sin(27rx)), x<0 , . 

1 = \ 300 x (1 + 0.1 sin(27rx)), x > 0 ' ^ 

The densities for each species can be solved by LTE condition. We add a small perturbation 
of the velocity over the whole domain 

u' = u + 0.05 x sin(7rx). (52) 

The computation stops at time t = 0.01. We remark that the solutions were computed on a 
larger domain [—6, 6] but truncated on [—5, 5] for not considering the effects by the boundary 
condition. Figs. 3 and 4 show the densities, temperatures, velocities and pressures by the 
balanced WENO-LF scheme and the regular WENO-LF scheme with 100 cells. The reference 
solution is computed by the WENO-Roe scheme with 1200 cells. From the Figs. 3 we can 
see that the balanced WENO-LF produces a more accurate result than the regular WENO- 
LF scheme. The regular WENO-LF scheme has a discrepancy from the reference solution on 
the waves and it cannot capture the small wave close to the shock. Unlike the density and 
temperature, the velocity and pressure are constant at the initial time. Thus it is more clear to 
see the difference between the balanced WENO-LF and the regular WENO-LF on the velocity 
and pressure results. From Figs. 4, we can see the results by the balanced WENO-LF are 
undistinguished from the reference solution. However, the regular WENO-LF produces large 
oscillations due to the truncation errors. 

Figs. 5 and 6 show the results by the central filter schemes. The considered central fil- 
ter schemes here are CEN2TVDfi and CEN4TVDfi with 200 cells, and CEN6WEN05fi and 
CEN8WEN07fi with 100 cells. All the well-balanced central schemes can capture the small 
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Figure 3: ID stationary contact discontinuity problem by WENO-LF schemes with 100 cells: 
left: density of O 2 ', right: temperature. 


perturbations very well. Figs. 7 and 8 show the results by PC filter schemes. The consid- 
ered PC filter schemes here are PC2TVDfi, PC4TVDfi, PC6WEN05fi and PC8WEN07h with 
200 cells. We remark that the PC scheme we implemented seems losing accuracy due to the 
boundary or the time discretization. Thus for some PC filter schemes we use a more refined 
mesh than the same order of the central filter schemes. But we will try to fix it in the future. 
The well-balanced PC filter schemes can capture the small perturbations well except a small 
wave around x = — 2 in the pressure plot (right subfigure of Figs. 8). This is because the 
well-balanced filter schemes are almost well-balanced everywhere except the transition areas 
between the filter part and non-filter part which are not well-balanced (see Sec. 5). The small 
wave will be resolved when the mesh is more refined. 


6.2.4 ID shock/turbulence interaction problem 

The fourth example is a ID Shu-Osher problem for reacting flows on the domain [—5, 5]. Initially 
a shock is located at x = —4. The shock is moving at the speed 500m/s to the right. The right 
state of the flow consists two parts, the first part is a constant equilibrium state from -4 to 1 
and the second part is an oscillatory equilibrium state from 1 to 5 with sine waves in densities 
and temperature. The conditions are given by 


{Tr,pr,vr) 


(500,20,0), a;e [-4,1] 

(500 x (1 + 0.1 sin(27ra:)), 20, 0), x G [1, 5] 


(53) 


Given the temperature and pressure, the densities for each species at the LTE state can be 
uniquely determined by the LTE condition. The left equilibrium state can be calculated ac- 
cording to the Rankine-Hugoniot jump condition. 

Since the right state of the flow is a zero velocity LTE state, the well-balanced schemes can 
solve it with machine round-off errors. If we add a small perturbation of the velocity all over 
the right state 

u! = u+ 10 -3 x sin(7ra;), x G [—4, 5], (54) 
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Figure 4: ID stationary contact discontinuity problem by WENO-LF schemes with 100 cells: 
left: velocity; right: pressure. 




Figure 5: ID stationary contact discontinuity problem by central filter schemes: left: density 
of O 2 ; right: temperature. 
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Figure 6: ID stationary contact discontinuity problem by central filter schemes: left: velocity; 
right: pressure. 




Figure 7: ID stationary contact discontinuity problem by PC filter schemes: left: density of 
O 2 ; right: temperature. 
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Figure 8: ID stationary contact discontinuity problem by PC filter schemes: left: velocity; 
right: pressure. 


the well-balanced schemes will be able to capture this small perturbation very well. 

At time equal to 0.012 seconds, the shock moves to x = 2. The shock first passed through 
the small perturbation of the velocity and then passed through one sine wave (at the region 

x e [1,2]). 

The reference solution is computed by WENO-Roe with 2000 cells. Figs. 9 and 10 show the 
comparison of the regular 5th order WENO-LF scheme and the balanced 5th order WENO-LF 
scheme on velocity and pressure with 100 cells. From the global views of velocity and pressure 
(the left subplots of Figs. 9 and 10), we hardly see any differences between these two schemes. 
However, if zooming in the region [—5,1], we can see the balanced WENO-LF scheme can 
capture the small perturbation very well whereas the regular WENO-LF scheme which is not 
well-balanced schemes cannot do that (the right subplots of Figs. 9 and 10). The density and 
temperature results by WENO-LF schemes are also shown in Figs. 11. 

The results by central and PC filter schemes are shown in Figs. 12, 13, 14, 15, 16 and 17. 
We remark that for the strong shock problem, the filter schemes cannot be essentially non- 
oscillatory around the shocks (although the oscillations are hardly seen from the global views). 
The well-balanced schemes are proposed to have advantages for the region close to the steady 
state, so we only focus on the regions away from the shocks (x £ [-5,-2]). The considered 
central filter schemes here are CEN2TVDfi and CEN4TVDfi with 300 cells, and CEN6WEN05fi 
and CEN8WEN07fi with 100 cells. The considered PC filter schemes here are PC2TVDfi and 
PC4TVDH with 300 cells, and PC6WEN05fi and PC8WEN07fi with 200 cells. Because for 
the moving shock problem, the filtered region is moving with the shock, and the switches 
between filtered and unfiltered will cause non well-balancedness, so the results of moving shock 
in this section are not so impressive as the results of stationary contact discontinuity in Sec. 
6.2.3. We use a more refined mesh for low order schemes CEN2TVDfi, CEN4TVDfi, PC2TVDh 
and PC4TVDfi. High order schemes CEN6WEN05h and CEN8WEN0711 with 100 cells have 
underresolved solutions in the crest and trough of the waves. PC6WEN05fi and PC8WEN07H 
with a more refined mesh 200 cells can resolve them better. 
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Figure 9: ID Shu-Osher problem of velocity results by WENO-LF schemes with 100 cells: left 
global; right: zoomed in. 
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Figure 10: ID Shu-Osher problem of pressure results by WENO-LF schemes with 100 cells 
left: global; right: zoomed in. 
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Figure 11: ID Shu-Osher problem by WENO-LF schemes with 100 cells: left: density; right: 
temperature. 



Figure 12: ID Shu-Osher problem of velocity results by central filter schemes: left: global; 
right: zoomed in. 
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Figure 13: ID Shu-Osher problem of pressure results by central filter schemes: left: global; 
right: zoomed in. 
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Figure 14: ID Shu-Osher problem by central filter schemes: left: density; right: temperature. 
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Figure 15: ID Shu-Osher problem of velocity results by PC filter schemes: left: global; right: 
zoomed in. 



Figure 16: ID Shu-Osher problem of pressure results by PC filter schemes: left: global; right: 
zoomed in. 
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Figure 17: ID Shu-Osher problem by PC filter schemes: left: density; right: temperature. 


6.2.5 A shock tube problem 

The last ID example is a shock tube problem. A diaphragm is located at x = 0 which separates 
the right chamber from the left chamber. The right chamber has cold air with pressure 0.6 x 
10 b N/m 2 and temperature 300K. The gas in the left chamber has high pressure 6 x 10 5 N/m 2 
and high temperature 3000K. Both gases are in LTE condition. The computational domain is 
[—5, 5] in the lab frame. 

The results are computed at t = 0.001. The solution is no longer in steady state. This 
example is to test the shock capturing ability of our well-balanced filter schemes. The numerical 
results of temperature, velocity and mass fraction of N 2 (from left to right) computed by 
CEN6WEN05fi, CEN8WEN07fi, PC6WEN05fi and PC8WEN07fi are plotted in Figs. 18, 
19, 20 and 21. As expected, the rarefaction wave, contact surface and shock appear in the 
temperature solution. Since the velocity is consistent through the contact surface, thus there 
are only rarefaction wave and shock appearing in the velocity solution. Furthermore, mass 
is conserved during the shock. No shock appears in the mass solution. All the considered 
well-balanced filter schemes can capture the shocks sharply without oscillations. 

6.3 Two dimensional numerical results 

As mentioned in the beginning, extending the well-balanced schemes to the zero velocity steady 
state of 2D reacting flow is trivial because the reacting term does not explicitly depend on the 
dimensions. In this section, similar well-balanced tests to 2D reacting flow will be performed. 

6.3.1 2D Well-balanced test 

Similar to ID, the first example is to check that our scheme maintains the 2D zero velocity 
steady state exactly. The 2D special stationary case 

T = 1000 x (1 + 0.2 sin(7r(x + y)))K, p = 10 5 N/m 2 , u = 0 m/s, (55) 
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Figure 18: Shock tube problem by CEN6WEN05fi: Left: temperature; Middle: velocity; Right: 
mass fraction of N 2 (CEN6WEN05fi with 300 points: dash-dot; Reference 1200 points: solid). 





Figure 19: Shock tube problem by CEN8WEN07fi: Left: temperature; Middle: velocity; Right: 
mass fraction of N 2 (CEN8WEN07fi with 300 points: dash-dot; Reference 1200 points: solid). 





Figure 20: Shock tube problem by PC6WEN05fi: Left: temperature; Middle: velocity; Right: 
mass fraction of N 2 (PC6WEN05fi with 300 points: dash-dot; Reference 1200 points: solid). 
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Figure 21: Shock tube problem by PC8WEN07fi: Left: temperature; Middle: velocity; Right: 
mass fraction of N% (PC8WEN07fi with 300 points: dash-dot; Reference 1200 points: solid). 


Table 2: L 1 relative errors for temperature by central/PC schemes at t = 0.01. 


N x N 

error 

error 

error 

error 


CEN2TVDfi 

CEN4TVDfi 

CEN6WEN05fi 

CEN8WEN07fi 

50 x 50 

4.14E-11 

4.14E-11 

4.21E-11 

4.30E-11 

100 x 100 

4.04E-11 

3.69E-11 

3.73E-11 

3.82E-11 


PC2TVDfi 

PC4TVDfi 

PC6WEN05fi 

PC8WEN07fi 

50 x 50 

4.13E-11 

4.20E-11 

4.24E-11 

4.29E-11 

100 x 100 

4.03E-11 

4.07E-11 

4.10E-11 

4.13E-11 


is considered. The computation is performed to t = 0.01 (about 2000 time steps for 100 x 100 
grid points) on the domain [—1, l] 2 . Table 2 shows the L 1 relative errors for the temperature T. 
We can clearly see that the L 1 relative errors are at the level of round-off errors, verifying the 
well-balancedness of the considered central filter and PC filter schemes for 2D reacting flow. 

6.3.2 2D small perturbation test 

The second example is again a small perturbation test but on a 2D steady state. The same 2D 
steady state solution Eq. (55) is considered. A small perturbation e = 10~ 3 x sin(7r(a: + y)) is 
added to the velocity in the x direction, i.e. , 

u = u + e. (56) 

The other quantities are kept unperturbed. The reference solution is computed by WENO-Roe 
scheme with 200 x 200 points. Fig. 22 show the contours of velocity by 2nd order central/PC 
TVD filter schemes at t = 0.01. The results by sixth and eighth central/PC WENO filter 
schemes are shown in Figs. 23 and 24. The ID cross-section results by PC filter schemes, central 
filter schemes and WENO-LF are shown in Fig. 25 left, middle and right subplots separately. 
We can see that our well-balanced filter schemes can capture the small perturbation in a coarse 
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Figure 22: 2D small perturbation of velocity results by filter schemes: e = 10 3 x sin(7r(a; + ?/)). 
Left: CEN2TVDfi 100 x 100 points; right: PC2TVDfi 100 x 100 points. 


mesh very well (especially for the schemes with order higher than 2 where only 40 x 40 is used). 
However, the WENO-LF produces large oscillations even in a mesh 100 x 100 (right subplots 
of Fig. 25). 

7 Concluding remarks 

In this paper the well-balanced approach is extended to the high order filter schemes in solving 
five species reacting flow in one and two space dimensions. Numerical examples are given to 
demonstrate the well-balanced property, accuracy, good capturing of the small perturbation 
of the steady state solutions, and the non-oscillatory shock resolution of the proposed well- 
balanced filter schemes. Because of the property of the zero velocity steady state solution of 
the reacting flow, the extension to any number of species and other reaction models is straight- 
forward. Future research will consider the non-zero velocity steady state and the advantages of 
well-balanced schemes to various steady state problems. 
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8 Appendix: Predictor-Corrector schemes and other spa- 
tial base schemes 

Samples of the high-order base schemes for F x can be of the following types. 

Central differencings: 
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Figure 23: 2D small perturbation of velocity results by central filter schemes: e = 10 3 x 
sin(7r(x + y)). Left: CEN6WEN05fi 40 x 40 points; right: CEN8WEN07fi 40 x 40 points. 



Figure 24: 2D small perturbation of velocity results by PC filter schemes: e = 10 3 x sin(7r(x + 
y)). Left: PC6WEN05fi 40 x 40 points; right: PC8WEN07fi 40 x 40 points. 
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Figure 25: Cross section of 2D velocity results at y = 0: e = 10 3 x sin(7r(x + y)). Left: central 
filter schemes; middle: PC filter schemes; right: WENO-LF. 


CEN4: 

F x ~ i2Ax ( Fj+ 2 ~ 8 ^ +1 + ^ F i~ l ~ F i - 2 )» ( 57 ) 

CEN6: 

Fx ” 6ok ( ^' +3 " 9 ^+ 2 + 45i ^ +1 " 45F ^- 1 + 9F i-2 - ^-3)- (58) 

Compact central differencings (Hirsh [8], Ciment and Leventhal [4], and Lele [12]). Here 


F x ~ (A x B x F)j , 

where for a fourth-order approximation 

(A x F)j = 4(^+1 + 4Fj + Fj-i), 

(B x F)j = 2 (Fj+i ~ Fj—i), 

and a sixth-order approximation 

( A x F)j = l(F) + i + 3 Fj + Fj- 1), 

{B x F)j = Aj(Fj + 2 + 28F}_|_i — 28Fj_i — F],_ 2 ). 

Predictor- corrector differencings: 

PC4: 

F^pF)' = 6A7 tf F j ~ & F j~ 1 + F j~ 2 ) i 

( — FFj + 8Fj+i — Fj+ 2 ) , 

PC6: 

F p F j = HaKx (37F} — 45F)_i 4- 9F)_2 — Fj_ 3 ) , 

D cFj = 30AJ (-37Fj + 45^+1 - 9^+2 -I- Fl^) , 

and PC8: 

F p F i = 42 UAX (533 Fj — 672Fj_i 4 - 168F)_2 — 32F)_3 + 3Fj—f) , 
D c Fj = 42qa^ ( — 533Fj + 672Fj + i — 168F)+2 4- 32Fj + 3 — 3F} + 4) , 


(59) 

(60) 
(61) 

(62) 

(63) 

(64) 
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where D p F is the PC differencing operator approximating F x at the first step (predictor step) 
and D C F is the PC differencing operator at the second step (corrector step). New forms of the 
upwind biased PC methods including compact formulations developed by Hixon and Turkel 
[9, 10] are also applicable as spatial base schemes. Interested readers should refer to their paper 
for the various upwind-biased PC formulae. The choice of the time integrators for these types 
of PC methods is more limited. For example, if second-order time accuracy is desired, then 
(62), (63) and (64) in conjunction with the appropriate second-order Runge-Kutta method 
are analogous to the familiar 2-4, 2-6 and 2-8 MacCormack schemes developed by Gottlieb and 
Turkel [5] and Bayliss et al. [1] . Here the first number refers to the order of accuracy for the time 
discretization and the second number refers to the order of accuracy or the spatial discretization. 
However, in this case one achieves the second-order time accuracy without dimensional splitting 
of the Strang type [25]. For higher than second-order time discretizations, only certain even 
stage Runge-Kutta methods are applicable. For compatible fourth-order Runge-Kutta time 
discretizations, see Hixon and Turkel for possible formulae. For example, the classical fourth- 
order Runge-Kutta is applicable provided one applies the predictor and the corrector step twice 
for the four stages, i.e., the predictor step for the first and third stages and the corrector step 
for the second and fourth stages. 

For the considered ID system with source term (2), the predictor-corrector scheme with 
2nd-order implicit explicit Runger-Kutta in time takes the form 

U (1) = U n - At.D p F{t n ,U n ) + AtS(t n ,U n ), (65) 

u n+1 = +U n )-AtD c F(i n+1 ,U w ) + AtS(t n+1 ,U n+1 ))/2, (66) 


The PC operators are modified at boundaries in a stable way by the sunnnation-by-part 
(SBP) operators [17, 16, 32], If D b is the standard pth order SBP for the centered difference 
operators, then the pth order PC operators are modified in this way 


F x 


at the predictor step and 


D P Fj > j — b + 1, ... ,N 
(■ D b -D c )Fj , j = l,...,b 


(67) 


f D c Fj , j = 1, . . . ,N — b 
\ (D b — D p )Fj, j = N — b+l,...,N 


(68) 


at the corrector step, where N is the number of grids and b is the number of points needed at 
boundaries. 
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